Flow equation approach to the linear response theory of superconductors 
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We apply the flow equation method for studying the current-current response function of electron 
systems with the pairing instability. To illustrate the specific scheme in which the flow equation 
procedure determines the two-particle Green's functions we reproduce the standard response kernel 
of the BCS superconductor. We next generalize this non-perturbative treatment considering the 
pairing field fluctuations. Our study indicates that the residual diamagnetic behavior detected above 
the transition temperature in the cuprate superconductors can originate from the noncondensed 
preformed pairs. 



I. INTRODUCTION 

Non-perturbative scheme of the flow equation method 
introduced by F. Wegner [l| and independently by K. 
Wilson with S. Glazek [2j proved to be a useful tool for in- 
vestigating a number of problems in the condensed mat- 
ter physics [|[ , mesoscopic systems [H , quantum optics [H 
and quantum chromodynamics [6|]. This procedure has 
been also recently applied to study the non-equilibrium 
transport phenomena of the correlated nanoscopic sys- 
tems Q- The main idea is rather simple and relies on 
a continuous process which, step by step, transforms 
the Hamiltonian to a diagonal or at least block-diagonal 
structure. One can use for this purpose various operators, 
depending on the subtleties of the discussed problem Q . 

Such continuous diagonalization scheme is reminiscent 
of the Renormalization Group (RG) technique [§[. They 
are similar with regard to the treatment of high/low en- 
ergy states (fast/slow modes). In initial steps of the con- 
tinuous transformation procedure mainly the most off- 
diagonal terms (i.e. high energy sector) are dealt with. 
Subsequently, the remaining parts closer to the diagonal 
are transformed. Since different energy scales are suc- 
cessively transformed/renormalized one by one the algo- 
rithm of flow equation method is relative to the family of 
RG formulations [sj . Let us remark that such techniques 
are in principle unrestricted by limitations of the usual 
perturbative methods. 

In this paper we: (1) formulate the current-current 
response function for the superconducting system using 
the flow equation method, and (2) extend such scheme 
to a state of the preformed pairs which above the criti- 
cal temperature T c loose the long-range coherence. Our 
study is motivated by the recent torque magnetome- 
try data of the Princeton group [l(| revealing the dia- 
magnetic features well preserved above T c in the lan- 
thanum and yttrium cuprate oxides. Similar indications 
have been also reported from the dc susceptibility mea- 
surements for Bi2.2Sri.8Ca2Cu30io+5 [HI and from the 
high-resolution SQUID data for Sm-based underdoped 
YBCO compounds [Hj]. Since the observed diamagnetic 
response is rather strong it can be hardly assigned to the 
Ginzburg-Landau fluctuations [l3j . 

Adopting argumentation discussed in the literature 



on the microscopic fl4l - fl9j and the phenomenological 
grounds [2(| [Hj] we consider the system consisting of the 
preformed local pairs (of arbitrary origin) coexisting and 
interacting with the itinerant electrons. Using the flow 
equation approach we analyze the diamagnetism within 
such scenario. Our study can be regarded as complemen- 
tary to the recent Quantum Monte Carlo (QMC) simu- 
lations for the same cooperon-fermion model [22j . It has 
also some resemblance to considerations of the supercon- 
ducting fluctuations beyond the gaussian approximation 
carried out in the t-J model [23j . 

We start by discussing the usual BCS model, treating 
it as a testing field for formulation of the linear response 
theory in terms of the flow equation method (readers 
less interested in the methodological details can skip this 
section). We next apply the same treatment to the mix- 
ture of electrons interacting through the Andreev scat- 
tering with the preformed local pairs. We determine the 
current-current response function and try to asses the 
diamagnetic response above T c . In summary, we point 
out the main conclusions and give a list of problems which 
might be of interest for further studies. 



II. BCS SUPERCONDUCTOR 



Let us first briefly illustrate how the flow equation pro- 
cedure determines the quasiparticle spectrum and the 
corresponding response function for the usual BCS model 
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describing electrons coupled to the pairing field Ak. We 
use here the standard notation for the creation (annihi- 
lation) operators c k (ck<r) and denote by £k = £k~ M the 
energies measured with respect to the chemical potential 
fx. Formally, Ak can be thought as the Bose-Einstein con- 
densate of the Cooper pairs Ak = — J^q ^k,q(c_ q j.c q -|-) 
which are formed by some attractive potential Vk,q < 0. 
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A. The continuous diagonalization 

The continuous diagonalization of the reduced BCS 
Hamiltonian (Ql has been considered in the original paper 
by F. Wegner [l[ and by several other authors [24], [25| . 
We briefly recollect the main steps of such procedure (see 
appendix A for the procedural details) which shall be use- 
ful for formulating the linear response theory. 

Following Wegner [l| we choose the generating opera- 
tor f]{l) defined by Eq. ([A3]) , which for the BCS model 
(fTJ) has the following structure 

fj(l) = 2£&(0 (A k (04t£i4 - K(l)c-H^t) • (2) 

k 

Transformation of the Hamiltonian H(l) proceeds as long 
as rj(l) is finite, which occurs until A k (Z) — > 0. This is 
achieved in the asymptotic limit I — > oo (jA4|. 

Substituting @ to the general flow equation (|A2[) for 
the Hamiltonian (TTJ) we obtain [25| 



dl 



dl 



K k (0 

lnA k (0 
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These equations pi4[) yield an exponential flo 



A k (l) = A k e" 



-4 f o ' rfi'[? k (r)] 2 



(3) 
(4) 

(5) 



and £k(0 = Cke 4 -^ d ' ' Ak ^ " , therefore the off-diagonal 
term A k (Z) vanishes in the limit I — > oo. Combining HI) 
we moreover notice that 4 {£ k (Z)) + |A k (/)| 2 } = which 
implies the invariance £ k (0 + |A k (/)| 2 = const. Due to 
A k (oo) = we conclude that the quasiparticle energies 
take the following BCS form 



^k = sgn(C k )^ + |A k |2, (6) 



where we introduced the shorthand notation for the 
asymptotic value £ k = lim;_ i . 00 £ k (Z). 



with the initial boundary conditions u k (0) = l, v k (0) = 0. 
Arranging the Z-dependent coefficients if front of c k ^ and 



c_ k ^ on both sides of the flow equation (|A7I) we find that 
du k (l) 



dl 
dl 



= -2^ k (l)At(0«k(0, (9) 
= 2&(l)A k (Z)u k (0- (10) 



From the equations (I9ll0j) we can see that the sum rule 
w k (Z) + v k (l) = 1 is properly conserved. To determine the 
needed asymptotic values we can rewrite © as = 

— 2£ k (Z)A k (Z)cZZ and substituting v k (l) = a/1 — u k (Z) we 
can analytically solve the integral J °° e?Z. In the asymp- 
totics we obtain the usual Bogoliubov-Valatin coefficients 



1 + & 



(11) 



Fourier transform of the single particle Green's function 



G a (k,iui) = /T 1 / dre'^ r G CT (k,r) 
Jo 



(12) 



(where j3 1 — ksT) takes hence the two-pole structure 



G CT (k, iw) 



,-,2 



iuj - 



(13) 



signaling the particle-hole mixing, characteristic for the 
BCS state. 



C. The linear response theory 

We now adopt the same procedure for studying a re- 
sponse of the BCS superconductor to a weak electromag- 
netic field A(r, t). In the linear response the induced cur- 
rent J(r, t) is assumed to be proportional to the pertur- 
bation, i.e. J(r, t) = -J dr' dt'K(r-~r', t-t')A(r', t'). 
Fourier transform of the kernel function consists of the 
diamagnetic and paramagnetic contributions [26j j 



B. The single particle excitations 

As an illustration how one can use this procedure to 
obtain various Green's functions let us derive the single 
particle excitation spectrum determined by G a (k, r) = 
— 3V(c kc r(T)c k(T ) , where T T denotes chronological ordering 
and r stands for the imaginary time. Applying to 
the flow equation (|A7[) for the creation and annihilation 
operators we infer the Bogoliubov ansatz 



c k f(0 = u k (l)c k1 - + w k (/)cL H , 
L (Z) = -« k (Z)c kt + u k (l)c> 



-kj.v 



kj. 



(7) 
(8) 



K a ^(q,uj) = — 5 a . p + e 2 U a . p (q,uj). (14) 

From now onwards, by a, /? we shall denote the Carte- 
sian coordinates x, y, and z. The paramagnetic term 
can be expressed by the (analytically continued) Fourier 
transform (|12p of the current-current Green's function 

n Q ^(q,r) = - (f T j qja (r)i_ q ,0>, (15) 

where the corresponding current operator j q = y + j_ 
consists of the spin f and J, contributions 



| e k,cr 5 k+q lCT 



(16) 



3 



with the velocity v k = ?i~ Vk£k- 

The standard way for computing the current-current 
response function (|15|) is based on the diagrammatic 
bubble-type contributions involving the particle and 
hole propagators and another contribution from the off- 
diagonal (in Nambu notation) single-particle propaga- 
tors. In this section we retrieve the standard BCS result 
[26j using the flow equation routine. 

To guess the relevant flow of the current operators we 



start by analyzing the initial (/ = 0) derivative 

= [ m wL (i7) 

\ / 2=0 

where jZ(l = 0) corresponds to the definition (TT51) . Using 
the generating operator ([2]) we find that 



I — ^— I = ±2 v k+f (&A k fi_ k ,_ a 8 k+qi<7 + £k+qA k+q c ko .ci. (k+q) 

V / 2=0 k 



(18) 



where the sign + (— ) refers to spin f (J,). Eq. (|18p unambiguously implies the following /-dependent parametrization 



jq(0 = Ek v k+a Mk,q(0c k , t a k+q.t + #k,q(0c-k,;cL (k+qU + 23 k, q (0c k ,f £ -( k +q)4 + -^(Oc-k^k+q.T 
jq(0 = Ek v k+f Mk,q(0c k ,4.Ck+q4 + B k,q(0C-k,tcL (k+q)>t - ^k,q(04,f 2 -(k+«lU ~ ^.qCOC-k^Ck+q^ 



(19) 



with the initial boundary conditions -4 k . q (0) = 1 and £>k, q (0) =l?k,q(0) =.Fk,q(0) = 0. All /-dependent coefficients have 
to be determined applying the ansatz (fl"9"|) in the flow equation (|A7|) for the current operators j q (/). On this basis we 
obtain the following set of equations 



<£A k , q (Z) 
all 

all 

rfPk,q(Q 

all 

dT^qjl) 
all 



-2 Kk + q(0A k +q(0^k,q(0 + £k (I) A k (1)^,(0] 

2 [£ k (Z)A k (Z)X> k , q (0 + ek+q(/)A k+q (/)J- k , q (/)] , 

2 [6c + q(0Ak+q(0A,q(0 - £ k (i) A k (Z)B k , q (/)] , 
2 Kk(0A k (0A,q(0 - a+q(0A k+ q(0Bk,q(0] • 



By inspecting (|20H23[) we can notice that 



d 



which implies 



J t [A,q(0 + S k ,q(0] = 2Kk+q(0Ak + q(/)-^(/)Ak(/)] [-F k , q (Z) - 23 k ,q(0] , 
| [2?k,q(0 - Jk,q(0] = -2 [a+q(0Ak + q(0 - £k(0M0] [A,q(0 + Sk, q (0] : 



| Hk, q (Z) + S k ,q(/)] 2 + | [2?k,q(0 ~ -F k ,q(0] 2 = 0. 



Taking into account the initial boundary conditions we hence obtain the following invariance 

[A,q(0 + £ k ,q(/)] 2 + [Pk,q(0 - ^k, q (0] 2 = 1 



(20) 
(21) 
(22) 
(23) 

(24) 
(25) 

(26) 
(27) 



valid for arbitrary /, including the limit I — > oo. Combining (|27p with the differential equations (|24l25p we exactly 
determine the asymptotic limit values 



•4 k . q + B kl , 



2> 



k,q 



k,q 



5 1 



A k + q Ak 



£k+ q 6c 



£ k +q£ k 
A k+q A k + £k+ q £k 



£k+q£k 



(28) 
(29) 
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where £k = \/£k + ^k- ^ n the same wa Y we also check that |/4.k,q(0 — #k, q (0] 2 — k,q(0 + ^qWf = 1 thereby the 
asymptotic values of all coefficients are found -4k. q = Uk£tk+q, £>k, q = Vk^k+q, 2? k,q = £kUk+ q , and -^k,q = UkVk+q- 

Since the transformed Hamiltonian H(po) is diagonal we can easily compute the current-current response function 
(IT51) expressing it through the particle-hole bubble diagrams (see the left h.s. panel in figure [T]). Finally it is given by 



n Q ,/j(q,w) = ^V k +f,aWk+^ 



A 



k,q 



B 



k.q 



f k,q — ^Tc,. 



1 - IFD (£k+q) - /i?D (£k 



/f.D (£k+q ) - /FA (£k ) 
1 



f 



1 



i^ + Ck+q-Ck «V-£ k +q+£k 
f 



W-£k+q-£k W + ?k+q + £k 



(30) 



where /fd(w) = [exp(w//cBT) + l] _1 is the Fermi-Dirac distribution function. We recognize that (|28l29p correspond 
to the usual BCS coherence factors (ttk+qUk + i>k+q^k) 2 and (wk+q^k — «k+ q Mk) 2 and thereby Eq. (l30|) rigorously 
reproduces the known BCS response function [26| . 



III. SUPERCONDUCTING FLUCTUATIONS 



A. Outline of the continuous diagonalization 



Numerous experimental data [HI [27l - [35j provided a 
rather clear evidence that the critical temperature T c in 
the underdoped cuprate superconductors (and similarly 
in the ultracold fermion gasses near the unitary limit 
t 38]) is not related to appearance of the fermion pairs 
but corresponds to the onset of their phase coherence. 
Upon approaching T c from above the short-range super- 
conducting correlations gradually emerge. For instance, 
the torque magnetometry [Toj and other measurements 
pll [l2j have detected the diamagnetic properties. 

To investigate the electrodynamic properties of the 
non-condensed preformed pairs we consider the model 

k,cr q 



k. P 



describing the itinerant electrons (c k ^ operators) coex- 
isting with the preformed pairs (bosonic Sq operators). 
They are mutually coupled through the charge exchange 
(Andreev-type) scattering. Such scenario (|3"Tj) has been 
considered by various authors in the context of high T c 
superconductivity (T^-tHI and for description of the reso- 
nant Feshbach interaction in the ultracold fermion atom 
gasses [36l - l38| . 

By Eq we denote the energy of preformed pairs mea- 
sured with respect to 2/z. Since in the superconducting 
state of cuprate materials the energy gap A k has a d- 
wave symmetry we furthermore impose the anisotropic 
boson- fermion coupling gk.p = .9 (cos k x — cos k y ). If one 
restricts only to the Bose-Einstein condensed pairs (i.e. to 
bosonic q = mode), then the model ([3Tj) becomes iden- 
tical with the reduced BCS Hamiltonian (TTJ) where Ak is 

related to the condensate — ^= 5k. -k- In what follows 
we shall consider an influence of the non-condensed pre- 
formed pairs on the current-current response function. 



Adopting again the Wegner's proposal [1| we choose 
the generating operator as fj(l) — [Hq(1), V B ~ F (I)], where 
Hq(1) stands for the free fermion and boson contributions 
whereas V B ~ F (l) denotes their interaction term. In ex- 
plicit form, such generating operator is given by 



fj(l) = 



-^=^a k , P (0 (& k +p< 

k,p 



■t =t 
kt C P ; 



h.C 



(32) 



where a k , P (0 = (£k(0 +£p(0 - £k+ P (0) 5k,p(0- Using 
([32]) in the flow equation for the Hamiltonian (|3I[) we 
obtain 39] 



dl 



]ng Kp (l) = -[&(l)+Ul)-E k+P (l)} 



(33) 



which implies an exponential diminishing of <?k,p(0 and 
guarantees its total disappearance in the asymptotic limit 
I — > oo. Simultaneously the fermion and boson energies 
are renormalized according to the flow equations (39j 



&(0 



$>k, q -k(0 5k, q -k(0 < (34) 



(35) 



where n^J B denote the fermion/boson occupancies of 



momentum k state. We have previously [39j, [40[ explored 
(analytically and numerically) the flow equations p3H35[) 
arriving at the following conclusions: 

a) the renormalized fermion dispersion £k develops ei- 
ther the true gap (below T c , when a finite fraction 
of the BE condensed bosons exists) or the pseudo- 
gap (for T C <T <T* , where T* marks the onset of 
superconducting- type correlations) , 

b) the long- wavelength limit of the effective boson dis- 
persion Eq is characterized by the Goldstone mode 
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(for T <T C ) whose remnants become overdamped 
in the pseudogap regime (above T c ), 

c) the single particle spectral function of fermions 
(see Appendix B) consists of the Bogolubov-type 
branches separated by the (pseudo)gap and these 
features remain preserved up to T* . 

More recently [l9[ we have also investigated evolution 
of the k-resolved pseudogap considering two-dimensional 
lattice dispersion with the nearest and next-nearest 
neighbor hopping integrals realistic for the cuprate super- 
conductors. For temperatures slightly above T c we have 
found that the pseudogap starts to close around the nodal 
points restoring the Fermi arcs, whereas in the antin- 
odal areas the pseudogap practically does not change. 
Upon a gradual increase of temperature the length of the 



Fermi arcs linearly increases, in agreement with the ex- 
perimental ARPES data [ijj. Similar conclusions have 
been achieved for the same model (|31|) from theoretical 
studies based the conserving diagrammatic approach [42j . 



B. Diamagnetism due to the preexisting pairs 

Following the guidelines discussed in section II. C we 
can now formulate the linear response theory for the 
model (|3ip . focusing on the role played by the non- 
condensed q^O preformed pairs. 

To impose the corresponding parametrization of the 
current operators j q we again start from the initial 
derivative (fT7|). Using the generating operator (|32|) we 
obtain 



dm 
di 



i=0 



where — (+) refers to the spin f (4)- Analyzing (|3"6"]l we deduce the following general structure of the Z-dependent 
current operators 

jq(0 = £ V k+§ (^k,q(0Ck,t ak +q.t + S k,q(0C-k,icl (k+q) i 
k 

+ ( P k,p,q(OSk+pc3 c ,tCp- q ,4. + -7 r k, P ,q(/)fo| c+p Cp^c k+q!t ) ) (37) 
p / 

with the initial values -4 k , q (0) = 1 and Bk,q(0) = £> k , P ,q(0) = J"k, P ,q(0) = 0. The operator j q (/) is given by expression 
analogous to (|37l) with 2?k,p,q(0 replaced by — J^.p q(0 & n d vice versa. Let us remark that taking into account only 
the BE condensed pairs & k + p = 6q ^p,-k we would come back to the ansatz (fT9|) reproducing the BCS solution. 
After a somewhat lengthy but rather straightforward algebra we derive the following set of the flow equations 

= £ [a k+ q, p _q(0© k ,p,q(0 «_ q + n£ +p ) + a k , p (0-Fk, P) q(0 « + < +p )] (38) 
p 

dBk ^ = - [a k:P (02 3 -P,-k, q (0 (rip + n k+p ) + afc+q.p-qCO-^-P ,-k,q(0 (n p _ q + ng +p )] (39) 



p 



rfPk,p,q(/) 
dl 

<^k,p, q (0 



-ak+q,p-q(/)^k,q(/) + a k , p (Z)£-p,q(0> ( 40 ) 
tjl -a kj p(0^k,q(0 + «k+q,p-q(0S-p,q(0. (41) 

For deriving the equations (|38I39[) we used the following approximations 

Kh' cJ,,V,« - ^k.k'^k c P , CT c p /^ + Kh' S p,p' n p - S k ,k' n k S p,p' u p ( 42 ) 
cp.fV.t c P)4 .c p ',4. - Kv n k c p ,4.c P ',4. + c p ,tV,t 5 P,P' n p - ^k.k'nk *P,P' n p (43) 

neglecting the higher order products SX SY of the fluctuations 5X = X — (X) , where the corresponding observables 
for the case of Eq. are defined by X — 6 k 6 k <, Y — c p o .c p / i<T and for (H31 by A = c p ^c p >^, Y — c p ^c p '.^. Such 

truncations 03]) enable us to satisfy the flow equation ^ j q (0 = ^(0>Jq(0 using the parametrization (|3"7) 
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k + q 



k-q 



q,+ k 



FIG. 1: (color online) Contributions to the current- current response function from the particle-hole bubble (the left h.s. panel) 
and from additional terms involving the finite momentum boson propagators (the middle and right h.s. panels). Vertices are 
expressed by the corresponding momentum components of the asymptotic values I oo for the coefficients used in Eq. (|37() . 



imposed on the current operators jq(0- Otherwise, if we introduced these neglected terms XY to the Z-dependent 
operator jq(0 they would induce even more complex structures arising from the commutator f}(l), jq(0 and formally 

the flow equation (|A7|) could never be obeyed (except for only the exactly solvable cases). The truncations (|42I43[) or 
similar, represent thus a necessary compromise in which the flow equation technique deals with the physical problems 
which are not exactly solvable [l|, ■ 

Finally, let us determine the current-current response function (|15|) . keeping in mind that the statistical averaging 
is feasible with respect to H(oo). Using the ansatz ([57)) we find the response function 



((jq,ai j— q,/3 



k,p \ a 

+ -4k,qSp i _ q ^((cj ii0 .Ck+q, ( 7; C_p j0 -C_(. p _ q ^ CT )) 
a 

+ ^k,q^p,-qy^((C-k,jcl( k+a ) „; C p cr Cp_q i0 -)) 

a 

+ S k , q 6p,_q^((c_ k , c7 C t _ (k+q)!cr ;C_p, (T C t _ (p _ q)!(T )) 

a 

- ^k,k',q^p,p',-q((6k+k'C k ,t a k'-q,;^P+p'Cp',iCp-q,t)) 
k',p' 

- ^ P ,p',q^k,k' -q<(&p + p/Cp',4.c p+q>t ;S p+p/ Cp/4Cp +q , t )) 
k', P ' 



(44) 



where G P , P ',q = I'p.p'.q — -^p.p'^q and we used the abbreviation ((Oi; 0%)) = —{T t Ox{t)C>2) h(oo)- These contributions 
(|44| arc depicted graphically in figure [1] Vertices denoted by the filled squares correspond to the asymptotic value Q 
whereas the filled circles represent A and/or B. 

Performing the Matsubara summation for the particle-hole convolutions (left panel in figure [lj and the double 
Matsubara summation for the diagrams involving one bosonic and two fermionic propagators we obtain the following 
Fourier transform of (l44l) 



-^k.a-^ 



k 

+ -4-k,-qZ?k,q + $k,q#k+q,-q 
+ ^ £k,-k',qgk+q,-(k'+q),-q 



k,q-^k+q, — q "t" --4k. q^— k, — q 



(45) 



fFD (£k+q) — JfD (£k 
1 ~ fFD (£k+q) - fFD (£k' ) 
l-/FZ)(£k'+q)-/FD(£k) 



1 



1 



W + ^k+q-^k W-Ck+q + Ck 

/Bg(^k-k') - /,Bg(£k+q + £k') 
W-(| k +q + |k' - -Ek-k') 

/Bg(-Ek-k') — /fiEKk'+q + £k) \ 
^+(fk'+q+|k - #k-k') J 



where /b_e(w) = [exp(ui/kBT) — 1] 1 is the Bose-Einstein distribution. The function (1431) in a straightforward manner 
generalizes the previous BCS form (|5Uf and is the central result of our study. 
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The d.c. diamagnetic properties of the system depend on the static value of the response function. In our present 
case it is given by 



n Q , Q (q,W = 0) = X! W k+2, Q S 2 A,q4+q ,-q + A:,qS-k,-q + A-k -q^k.q + Sk.q^lc+q ,-q 



Sfd (Ck+q) - ]fd (£k) 



£k,-k',q£k+q,-(k'+q),-q /flB^k-k') ~ /sE^k+q + Ik' 



fBE(Ek-k') — /s£;(£k'+q + £k) 
I 



Ck+q-Ck 

l-/FD(Ck+q)-/Fp(^ k ) 
£k-k' - (6k+q + £k') 

-Bk-k' - (£k'+q+£k) / I 



(46) 



For temperatures below T c (when a finite fraction of the 
Bose-Einstein condensed pairs exists) the main contri- 
bution in the expression (|46|) comes from k' = k terms. 
Under such conditions Eq. (|46|) becomes identical with 
the BCS solution consisting of: a) the superfluid fraction 
[i.e. the term in Eq. (|5P|) proportional to the coherence 
factor (uk+q^k — £'k+qU k ) 2 ], and the other b) normal con- 
tribution from the thermally excited quasiparticles, i.e. 
the term proportional to (M k + q u k + '5 k +qi> k ) 2 [E3|. 

Above T c (but fairly below T* ) a considerable amount 
of the preformed pairs occupies the low momenta states 
Eq^o therefore expression (14TH) becomes reminiscent of 
the above mentioned BCS components in the response 
function. 



IV. ITERATIVE SOLUTION 

To explore the physical aspects related to the current- 
current response function (|45|) we adopt an iterative 



method for solving the coupled flow equations (|38M1[) . 
Such scheme allows for an approximate estimation of 
the introduced /-dependent parameters. Following [39| 
we make use of the fact that the dominant renormal- 
ization affects the boson- fermion coupling gk, P (l) which 
ultimately vanishes in the asymptotic limit I — > oo. Ne- 
glecting the simultaneous renormalization of the fermion 
£ k (£) ~ £ k and boson energies Eq(l) 
following solution of the flow equation ([33 



Eq we obtain the 



e -(£k+£p--E k+p )' I 



(47) 



Substituting this result ([47"1) to the flow equations (|34I35|) 
we might in turn update the energies and the routine can 
be continued, at each iterative level providing a better 
and better estimation for the renormalized quantities. 

In this section we apply such scheme to the flow equa- 
tions (|38M41[) . restricting ourselves to the lowest order 
solutions based on Eq. (|47|) . We start with the initial 
values of the coefficients -4 kjq (/) ~ -4 k , q (0) = 1 and 
Bk,q(l) — Sk,q(0) = substituting them in the right h.s. 
of the flow equations (|40I41|) . Using Eq. (g7J| we analyt- 



ically solve the simplified equations (|40I41|) obtaining 



Pk,P,q(0 



■^k,p,q(0 



.9k+q,p— q 



3 -(4k+q-Kp-q--Ek+p) 2 i _ I 



5k,p 



Ck+q + Cp-q — ^k+p 

f,— (?k+Cp — -Ek+p) I \ 



ik + £ P - E k + P 

Their asymptotic values are given by 

9k+q,p-q 



^k.p.q 
*^"k,p,q 



£k+q + Cp-q - Ek^ 
ffk,p 

Ck + Cp - Ek+p 



(48) 
(49) 

(50) 
(51) 



Using the /-dependent coefficients £>k, p , q (0 and J"k, p ,q(0 
we can next determine Ak,p(l) and £> klP (/)- Substituting 
(I48I49[) in the right h.s. of (138139)) we obtain the following 
asymptotic values 



A 



k.q 



2 ^ 
p 



p-q 



\ n p +«k+ P ) Iflk-P 



(£k + Cp - ^k+p) 
k+p) l5k+q,p-q| 2 



and 



Sk,q 



n p + "k+p 



(Ck+q + £p-q - Ek+p) 



^2 -9k,p.9k+q,p-q 
P 

B 



(52) 



(53) 



1 



^k,i 



^k+ 



p-q 



q,p q 



l k+p 



Xk,\ 



y2 

^k+q,p-q / 
A k +q !p _c 



Xk, 



k+q,p-q 



X lp 



^k+q.p-q, 



where Xk, P = S,k+^ P - Ek+ p . 

Since eventual diamagnctism is determined by the long 
wavelength limit of the static response function (|46|) we 
focus on q = values of the coefficients. Examining the 
q+0 limit of the asymptotic values (|50I51[) we notice 
that the superfluid vertices vanish 



G 



k,p,q : 



: -Dk,p,q — ^k,p,q — > 



(54) 
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-1 q x a l 

FIG. 2: (color online) The renormalized boson energy Eq 
obtained for three representative temperature regions: T > 
T* (top panel), T* > T > T c (middle plot) and T c > T 
(bottom panel). The dashed lines show the level 2fi(T). 

and (similarly to the BCS treatment [5(|) we are left only 
with the normal component of the paramagnetic term 

hm n Q ,(q,0)-2V^[4o+6k,o] 2 — °- (6c) ■ (55) 

At high temperatures (in a normal state) the disper- 
sion £k is nearly identical with the bare energy £k — fi 
therefore Eq. (1551) cancels out the diamagnetic term of 
the response kernel K at p(q —> 0, 0) and consequently the 
system does not show any diamagnetic features. On the 
other hand, in the superconducting state the single par- 
ticle excitations £k are gaped and at low temperatures 
-J-/Fr>(£k) ~ ~ <5(6c) therefore the paramagnetic contri- 

"sk r 1 

bution (|55|) vanishes 50] . One then obtains a perfect dia- 
magnetism with the characteristic London penetration 
depth XJ 2 = ne 2 /m. Between these extreme regimes we 
can expect an intermediate behavior. In particular, for 
temperatures T c < T < T* the single particle fermion 
spectrum becomes partly depleted around the Fermi en- 
ergy so the paramagnetic term (|55[) would no longer be 
able to compensate completely the diamagnetic contri- 
bution generating a fragile diamagnetism. 

For some quantitative illustration of this behavior we 
have analyzed temperature dependence of the superfluid 
fraction n s (T) defined by the relation (2(| 

e 2 n (T) 

■4(q^0,0) = - sV ^ (q^O.O). (56) 

TO 

In order to determine n s (T) we substituted the para- 
magnetic term (f55"j) to the kernel function K a ,p(q — > 
0,0) and applied the coefficients (|52l53p . simplifying the 
fermion and boson concentrations by n£ ~ /fD(?k) and 




0.01 1 1 ■ 1 1 

0.05 0.1 0.15 0.2 

k B T/D 

FIG. 3: (color online) Temperature dependence of the effec- 
tive boson mass m B obtained for the initially discrete energy 
level Eq(l — 0) = const. Our results resemble the QMC data 
shown in Fig. 9a of Ref. J^. 

n q ~ fBE{Eq). Furthermore, we replaced all ener- 
gies by the renormalized values £k and Eq to account 
for the iterative feedback effects. Following the previ- 
ous study [l^| we have selfconsistently determined these 
renormalized energies £k> Eq by solving numerically the 
flow equations <|34I35|) for the tight-binding lattice model 
ck = — 2i [cos ak x + cos ak y ] — 2t z cos ck z assuming re- 
duced mobility along z-axis t z — O.lt. Initially (at ^ = 0) 
we have assumed bosons to be localized. To establish 
some correspondence with the recent QMC studies [1H 
we have used the same total concentration of carriers 0.16 
and imposed the coupling g — 0.2D (where D — 8t). 

The important changeover of the boson dispersion Eq 
upon varying temperature is shown in figure [2] We no- 
ticed that below some characteristic temperature ksT* ~ 
0.05D there occurred a considerable reduction of the in- 
plane boson mass, defined as d 2 Eq/dq 2 = h 2 /m B . Its 
temperature dependence [compared to the bare planar 
mass of fermions ~ Ti 2 /2ta 2 ] is illustrated in figure 
[3] The mentioned suppression of the boson mass be- 
low T* coincided with appearance of the pseudogap in 
the fermion spectrum near /j, - this property has been 
discussed at length in our previous studies [H, H(| where 
we formulated the flow equation procedure for the present 
model ([3T]) . Below the other temperature ksT c ~ 0.026-D 
the Bose-Einstein (BE) condensate appeared in the sys- 
tem and simultaneously the parabolic dispersion evolved 
into the collective sound-wave mode Eq oc |q| (see the 
bottom panel in Fig. [5]). 

Evolution of the effective boson and fermion spectra re- 
vealed a substantial influence on the superfluid fraction. 
In figure [?] we show the temperature dependence of such 
n s (T). Below the temperature T* (for here chosen set of 
the model parameters T* ~ 2T C ) we observed a gradual 
buildup of the superfluid fraction. Passing below T c the 
superfluid fraction exhibited a further, stronger enhance- 
ment manifesting an onset of the long-range phase coher- 
ence caused by appearance of the Bose Einstein conden- 
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0.02 0.04 0.06 0.08 0.1 
k R T/D 



FIG. 4: (color online) Temperature dependence of the su- 
perfluid fraction n s (T) normalized to the total fermion con- 
centration n. We can notice a broad regime where the fragile 
Meissner effect is caused by the superconducting fluctuations. 



sate of pairs. At still lower temperatures, i.e. deep in the 
superconducting state T <C T Cl we observed some flat- 
tening of the superfluid density, rather than the expected 
linear dependence n s (T — »• 0) ~ n s (0) — aT typical for d- 
wave superconducting systems with the Dirac-like excita- 
tions around the nodal points [5L, 52]. Presumably these 
low temperature results indicate that we are not correctly 
evaluating the transverse Fermi velocity va and/or the 
longitudinal one v/ therefore the proper linearity coef- 
ficient a = [2ln(2)/n]vf /va [EH is missing. We also 
suspect that in our computations we might apparently 
overestimate the role of antinodal areas, where the ma- 
jority of bosons is effectively gathered for T — > (see 
figure 4 in the Ref. 19]). This artificial low-temperature 
dependence of n s (T) needs a more careful investigation. 

Summarizing this section, we have obtained the super- 
fluid density n s (T) which clearly indicates a fairly broad 
temperature regime T c < T < T* of the Meissner rigid- 
ity appearing due to the superconducting fluctuations. 
Such fragile diamagnetism originates solely from the non- 
condensed preformed pairs, as has been previously sug- 
gested by several authors (43l - [45j . We have evaluated 
n s (T) by means of the new nonperturbative method and 
determined the current-current response function (|45[) 
solving iteratively the set of flow equations. Our results 
seem to qualitatively capture the experimental data of 
N.P. Ong group ond other recent measurements [53$ . 



such treatment onto the mixture of the non-condensed 
(preformed) pairs interacting with the mobile electrons 
through the charge-exchange Andreev scattering. Wc 
have determined the contributions (see Fig. [T|) to the re- 
sponse function (|45|) . where the vertices are expressed by 
the corresponding flow equations (|38M4ip . 

The central result ([43]) of our study generalizes the 
BCS current-current response function [261 ] taking into 
account the residual diamagnetic effects originating from 
the finite-momentum preformed pairs. Such effects are 
studied here in an alternative way than the perturba- 
tive corrections due to the Aslamasov-Larkin and Maki- 
Thompson diagrams [l3[ . In our approach the fluctations 
enter the current-current response function through the 
convolution of one boson and two fermion propagators 
(see figure 1) instead of the higher order convolutions typ- 
ical for the standard diagrammatic study. In the present 
formalism the influence of fluctuations affects the vertex 
functions which have to be determined from the asymp- 
totic solution of the flow equations (|38H41[) . 

In the static, long wavelength limit we find a clear 
evidence for the pronounced diamagnetic contribution, 
which might be relevant to the experimental data ob- 
tained for the underdoped cuprate materials in the lower 
part of the pseudogap state [10l4l2| . Our study is consis- 
tent with the recent Quantum Monte Carlo results for the 
same model [22j ■ In both cases the residual diamagnetism 
originates from the preformed pairs whose mobility con- 
siderably increases below T*. Similar ideas concerning 
the non-condensed pairs have been emphasized by sev- 
eral other authors [23l. [43l448j . 

In order to see a more specific relation of this treat- 
ment to the cuprate oxides one should solve numerically 
the set of flow equations (|3"81 -01]) for the realistic model 
with nearest neighbor and next nearest neighbor hopping 
integrals. Another issue (not addressed here) concerns 
the doping effects which should affect the Fermi surface 
topology and influence the populations of the fermions 
and the preformed pairs [22j, |39( . It would be also worth- 
while to solve the flow equations (|38M4ip fully selfconsis- 
tently and investigate the electrodynamic properties us- 
ing the response function K a ^(q, uj), which generalizes 
the standard BCS result. 

We hope that the present formulation of the linear 
response theory by means of the flow equation method 
could stimulate further studies of the many-body effects 
in various subdisciplines of physics. 



V. SUMMARY 

We have addressed the linear response of the electron 
system with the pairing instabilities using nonperturba- 
tive framework of the continuous unitary transformation 
technique For the case of the Bose-Einstein con- 
densed pairs we have analytically derived the BCS result 
(|3H)l . which in the static and long wavelength limit ac- 
counts for the Meissner effect. We have next extended 
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Appendix A: Methodology of the flow equations 



Appendix B: Effective quasiparticles above T c 



We give here a brief outline of the continuous canonical 
transformation for arbitrary Hamiltonian of the following 
structure 



H = H + H 1 



(Al) 



where Hq denotes the diagonal part (for instance it can 
be the kinetic energy of particles) and Hi stands for 
the off-diagonal term (i.e. interactions or any perturba- 
tions) . Upon continuously transforming the Hamiltonian 
H(l) — W (l)HU(l) the ^-dependence (flow) is governed 
according to the following differential equation PJ 



dH(l) 
dl 



W),H(l)] 



(A2) 



dmiu-i(i)_ 



where the generating operator #)(/) = —jt- 

Choice of 77 (Z ) is usually dictated by the specific physi- 
cal situation. One of the possibilities, suggested by Weg- 
ner [lj, is 



f}(l)= J3b(Z),ffi(J) 



It has been proved that (|A3|) guarantees 
lim Hi(l) = 



(A3) 



(A4) 



provided that no degeneracies are encountered. Several 
alternative proposals for rj(l) capable to deal with the 
divergences [2| , the degenerate states or various other 
advantages has been discussed in the monograph Q. 
To carry out the statistical averages of the observables 



(d>£ = Tr{e-^d}/Tr{e-^} 



(A5) 



(where j3 1 = ksT) it is convenient to use the invariance 
of trace on the unitary transformations 



Tr 



{e-^6}=Tr{e-^«o W } ) ( A6 ) 



where 6(1) = W(Z)6W _1 (0- The '-dependence of 6(1) is 
imposed through the flow equation [l| 



dO(l) 
dl 



= W)M)\ 



(A7) 



similar to (|A2j) for H(l). Since the Hamiltonian H(l) be- 
comes diagonal for I — > 00 the easiest way to compute 
the trace (|A6[) is with respect to H(oo). This however 
requires, that simultaneously with the continuous diago- 
nalization of the Hamiltonian one has to analyze the flow 
of other physical observables 6 — > 6(1) — > 6(00). 



To determine the single particle excitation spectrum 
for the model (I3T1) we have to transform the individual 



operators cj^ (l) = U (l)cj!l U 1 (l) which is a bit tricky be- 
cause U(l) is not known explicitly. Following the scheme 
outlined in section II. B and using the operator fj(l) cho- 
sen in the form (J32j) we deduce the following ansatz for 
the fermion operators [40j 



Ckt(0 
1 



Wk(0 Ck-f + Wk(0 c^ ki 

[wk,q(0 ^ q C q+ kt + Uk,q(0 &qC q - k4 . 



-kj. 



(l) = -v* k (l)c kt + K(l) cl 



(Bl) 



(B2) 



1 



A^ 

q^O 



[-«k, q (0 K 



q c q+kt T U k q 



(0 M 



where itk(0) = 1 and the other coefficients are vanishing 
at 1 = 0. These Z-dependent coefficients can be derived 
from the equation (IA7|) for c^(l) and cL k ^(l) operators. 
The corresponding set of flow equations reads (lol ] 



du^l) 
dl 



dvk(l) 
dl 



du k ^(l) 

dl 

dVk,q(l) 
dl 



N 

q#0 



-k,k 



*q-k,k 



(I) V k (l) 



(I) 



ak,k(0 u k (l) 
a k, q +k(0 (n% + 

q#0 

"q-k,k(0 Uk(0) 

- a k , q +k(0 u k(0- 



(B3) 

-kj Wk,q(0, 

(B4) 

-kf) U k,q(0, 

(B5) 
(B6) 



They are additionally coupled to the flow equations 
[33]) because of the terms ak,k'(0- If one neglects the fi- 
nite momentum boson states [when Mk, q (0 = = i>k,q(0] 
these equations can be solved analytically [4(3] , reproduc- 
ing the standard BCS case discussed in section II. B. The 
case of q ^ bosons is more cumbersome. We have 
previously studied such problem numerically (3tl l40j . 
in particular considering also the 2-dimensional square 
lattice [lH with the tight-binding dispersion £ k (0) = 
— 2t [cos(fc x a) + cos(k y a)] — At' cos(k x a)cos(k y a)— \i assum- 
ing the initial discrete energy £q(0) =Eq and fixing the 

total charge concentration 2 J2q n q + J2k ( n kt + n k4.) ■ 
One of the valuable results obtained from such for- 
malism concerned the pseudogap regime. Since above 
T > T c the condensate fraction is absent we can notice 
that (IB4IIB5P imply Uk(/) = = Uk,q(0- ^ n other words 
the parametrization (|B1[ IB2|) simplify then to 



Ckt(0 = Uk(0 Ckt + -^^"^(0 MUi ( B7 ) 

q^O 



1 
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-14 (0 = W k(0 ~ -J= Yl U k,q(0 ^qCq+kt, (B8) 



with the coefficients obeying 



— = ^ 2^ a q-k,k(0 + ™q-kj «k,q(<). 

^k,q(Q 



q#0 

ak, q +k(0 w k(0- 



Under such circumstances thus find that the single par- 
ticle spectral function A(k,w) = — 7r _1 ImG (T (k, ui + i0 + ) 
takes the following structure 



i4(k,w) = \u k \ 2 S (w-|k) 



(B9) 



1 



-Jr S K 

q/O 



^q-H) |6k,q| 2 5(w + f q -k--Eq). 



The first part of (|B9|) describes the quasiparticle whose 
renormalized dispersion £k in the temperature regime 
T < T* c is discontinuous at the chemical potential (sig- 
naling a depletion of the low energy states). The second 
part of Eq. (IB9[) contributes the shadow to the above 
mentioned quasiparticle. These contributions constitute 
the characteristic Bogoliubov-type excitation spectrum 
which has been indeed observed experimentally in the 
yttrium [3l[ and lanthanum 32] compounds. 
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